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We are studying Runge-Kutta methods along complex paths of integration 
from a geometric point of view. Thereby we derive special complex time grids, 
which applied to the problem of integrating a linear autonomous system of or- 
dinary differential equations, can be used to achieve a classical super convergence 
effect. The approach is also adapted for arbitrary ODEs. Furthermore we draw 
a connection from our geometric reasoning to the class of composition methods 
with complex coefficients. Thereby, our main goal is to introduce a new point of 
view on these methods. 
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' 1 Introduction 

O ' 

Numerical integration of ODEs is in essence a well understood subject and there is a wide 
range of methods that are applicable for different flavors of the problem. This paper aims at 
adding some novel twist to the subject: superconvergence in the presence of certain complex 
time grids and a new viewpoint on composition methods with complex coefficients. Usually 
numerical integration of ODEs is embedded into a context where the idealized time variable 
is considered to be a continuously flowing real entity. Numerical integration methods (like 
for instance variants of the Runge-Kutta approach) approximate this continuous flow by 
time steps of finite resolution. Usually the finer the grid the smaller the global approximation 
error. Adaptive step width algorithms aim to produce a compromise of numerical accuracy 
and computational effort. Near numerically instable situations, adaptive stepwidth algorithms 
tend to introduce many additional grid points. In fact, in a scenario of time flow within the 
real numbers there is no way to bypass these numerically critical situations. For instance, 
imagine a version of the restricted two-body problem where a planet moves around a fixed star 
under Newtons law of gravity. If the initial velocity of the planet points directly towards the 
sun, the system will necessarily run into a singular situation. Scenarios nearby this situation 
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will cause numerically difficult situations. In a naive adaptive step width approach usually 
many additional time steps are introduced when the planet swings closely around the sun. A 
setup that allows also complex time, leaves the possibility to circumvent the singular situation 
by introducing a time flow that makes a "detour" through complex values. 

Our research was initially motivated by the question of applying complex detours to resolve 
singularities in numerical integration for ODEs. The philosophy behind this approach has been 
successfully applied to other situations. For instance in the field of Dynamic Geometry (see for 
instance [H]) the second author successfully applied complex detours to avoid discontinuous 
behavior, which was known to be a notorious problem in this field |10j . In relation to making 
real-time numerical simulation available to the dynamic geometry program Cinderella, it 
was natural to study a complex detour approach, as well. 

Investigating complex detours for ODEs raises several interesting research questions. Be- 
sides many interesting implementation-related issues, in particular problems arise, which are 
related to the analytic character of the input and the output of the solver. 

• Analytic continuation of the right side: Functions in a complex scenario require 
the proper treatment of different branches (for instance for y / TT7 and log(. . .)). To 
avoid path jumping (as mentioned for dynamic geometry) all function evaluations have 
to mimic an analytic behavior. Thus in the evaluations of intermediate results the 
"correct" branches have to be chosen. This requires the implementation of complex 
tracing strategies for the evaluation of the right side of the ODE if the derivative in the 
ODE has several branches. 

• Complex manifold of the solution: Not only the right side of the ODE may exhibit 
monodromy behavior of an analytic function. Also the solution itself may depend on 
the actual path chosen for the time variable. Even if the right side of the ODE is 
unbranched, it may happen that the solution at a concrete moment of time depends 
on the integration path, connecting the start time and the end time. Such effects have 
for instance been studied in [3], where also a connection of chaotic behavior and these 
monodromy effects is drawn. 

Perhaps due to these reasons (the additional implementation problems and the path depen- 
dence of the solution) in the literature there are only very few cases in which the paradigm of 
real time is broken [2]. Still we are convinced that the complex setup is worth to be studied 
both from a practical standpoint that asks for good approximations with small computational 
effort, as well as from a structural purely mathematical standpoint. In the first part of this 
article we will focus on a specific situation, which in a surprising way connects these two 
aspects. We will study linear ODEs, which are interesting since they form a perhaps simplest 
possible scenario, in which one can apply complex detours. For a linear right side one has just 
one single branch, so we can neglect the analytic continuation issues. Surprisingly, in this case 
carefully chosen complex time grids can help to reduce the global integration error by at least 
one order. Nevertheless, there is also a disadvantage of this approach. Roughly speaking, the 
effort of complex computations is six times higher than in the real-valued setup. Since there 
exist integration methods of arbitrary order, it is natural to ask the following question: 

Why should one study numerical integration along complex time grids to increase the order 
of a one-step method by one, if one could use a real-valued one-step method, which has already 
got the increased order? 
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The answer is, that the authors have not been faced with this question at the beginning of 
their studies concerning techniques to detour singularities, but for the reader, which is familiar 
with the field of numerical integration, it should be interesting that the naive and geometric 
approach of the authors has led to a new derivation of the class of composition methods with 
complex coefficients, apart from the usual way of solving some suitable set of order conditions 
over the reals or the complex numbers like done in [8l[9j[l3]. In addition to that, the authors 
also discovered some results concerning the convergence of Runge-Kutta methods. 

The paper is organized as follows. In Section 2 we present a simplest possible scenario 
where the super convergence effects studied in this article arise. This example (namely complex 
detours for x = x) will serve as a motivating paradigm for our further considerations. Section 
3 introduces the necessary setup of complex time grids in relation to Runge-Kutta methods. 
Section 4 is the main technical part of this article. We first deal with the problem under which 
condition a Runge-Kutta method applied to a complex path yields a real terminal point 
(Theorem 14. ip . After this we proof a main result of this article: One can increase the order 
of a convergence of a Runge-Kutta method applied to a linear ODE by choosing a suitable 
complex path (Corollary 14. ip . We gain a lot of geometrical insight in the structure of the path 
and can derive explicit criteria that have to be satisfied in order to obtain super convergence. 
These criteria are closely related to the multiplicative structure of roots of unity in the complex 
plane. Section 5 draws the connection of the complex detours to composition methods, which 
enables us to extend our method (at least in a certain sense) also to the case of a non-linear 
right-hand side. The condition equation for super convergence derived in Section 4 are closely 
related to the classical criteria for increasing the order of composition methods. A composition 
method is obtained by applying a integration method (say a Runge-Kutta integration) 
consecutively in a controlled way. The sequence of these applications can again be considered 
as one step of a more complicated integration method. Under certain circumstances (the 
above mentioned criteria) one can increase the order of the original method by at least one. 
Our geometric approach allows to interpret some of these composition methods entirely on the 
level of an underlying integration grid. Thus these composition methods simply correspond 
to a suitably chosen complex detour. We can even iterate this process and by this obtain 
(Section 5.3) composition methods of arbitrary high degree (also known as the "Yoshida 
trick"). Surprisingly, the corresponding paths that encode the composition structure exhibit 
a fractal structure. Finally at the end of Section 5, we illustrate our methods by the more 
sophisticated problem of computing the Arenstorf orbit. 

2 A motivating example 

In this section we would like to make the reader familiar with the subject of numerical inte- 
gration along complex paths (and its benefits) in an informal but hopefully self-explanatory 
way. In order to concentrate on the fundamental ideas, this and the follwing section is written 
in a very easy and self-containing manner. 

We choose the simplest possible initial value problem 



where / : R 2 — > R, (x,t) x, and (to,Xo) := (0,1). The corresponding analytic solution 
curve 99 : R — ► R, 1 1— > (p(t), is given by the well-known exponential function <p(t) := e l , Vt S R. 



x(t) = f(t,x(t)), x(t )=x 
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In order to compute <p(l), a standard approach is given by using an explicit Runge-Kutta 
method (eRKM) along an equidistant decomposition of the interval [0, 1]. 

In contrast to the real interval [0, 1], we now study the use of complex paths connecting 
and 1. As (H|) is a well-behaved3 initial value problem given by an autonomous linear first order 
differential equation with constant coefficients, t S 1 seems to be an unnecessary restriction. 
To be more precisely we note the following. 

Our initial value problem (pQ) could be written equivalently as the integral equation 



where 7 : [to, t] — > C, s 1— ► s, is a (^-curve. As the solution 93 has to be a holomorphic function^, 
Cauchy's well-known integral formula states, that ip(t) is independent of the detailed choice 
of 7. Every C^-curve 7, with starting point to and ending point t, yields the same result for 



With this "analytic" picture in mind, we have a look at the behavior of the explicit Euler 
method along the following complex time grid. We now take 



as a discretization of the upper complex half circle from to 1. Starting at the initial value 
xq = 1, the explicit Euler method generates n, not necessary real-valued, approximation 
values 



In the following, we compare the explicit Euler method along the mentioned complex time 
grid and the equidistant composition of the interval [0, 1], given by t j := ^, Vj G {0, . . . , n}. 

Figure Q] illustrates this construction in the real and complex case. Note that Figure [1] 
represents only the image of the corresponding Euler polygons (to-plane). Thus the real- 
valued construction degenerates onto the real line. Looking at the numerical results, 



For existence and uniqueness of a complex solution curve have a look at [12] . 

2 In the case of problem {TJ, there exists a unique entire function as solution curve for every choice of (to,xo) £ 
C 2 as initial condition. 



x(t) = xq + x(s)ds = xq + / x(z)dz, 




<p(t). 



_( e **(i-4)+l), VjG{0,...,n} 



Xj ss <p(tj) = e tj , Vj € {0, . . . , n}. 
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Figure 1: Real and complex approximation for n = 10. Thereby + and x represent the 
corresponding exact values of the solution cp along the upper complex half circle 
respectively the real interval from to 1. □ and o are the corresponding values of 
the explicit Euler method. 
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we point out the following observations: 

1. Both terminal points seem to be an approximation of e ~ 2.718281828 (as expected). 

2. The imaginary part of the complex terminal point seems to be zero. 

3. The complex construction yields more correct digits! To be more precise, the complex 
terminal point xio is more than 16 times closer to e than the real terminal point xio- 

Especially the last observation seems a worthwhile phenomenon to go into a deeper study 
of problem ((TJ) and related ones. In this article we will develop a theorerrH, which explains 
this phenomenon in the more general context of numerical integration of linear initial value 
problems with constant coefficients by the use of explicit Runge-Kutta methods (eRKM)s 
along complex paths. 

3 Compare Theorem 14.21 and Corollary 14.11 
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3 Complex flows and one-step methods 

Before we confine our attention to the encountered observations above, we have to fix the 
notation in this section. 

Definition 3.1 (CIVP). Let d G N, Q f C C x C d and / : 0/ -> C d , (t,x) >-> f(t,x), be a 
continuous function. Then 

x{t) = f[t,x{t)), (2) 

is called an explicit first order differential equation. A solution curve of ([2]) is a complex 
differentiable function ip : U — > C , i i— > with 

1. := { (t, y>(*)) |t G 17} C 0/ and 

2. <p(t) = f(t,<p(t)), Vtef/, 

where [/ C C is an open set. ([2]) together with a point (to,^o) £ ^/j the initial condition, 
is called a complex initial value problem (CIVP). A solution to a (CIVP) is a solution curve 
99 : 17 -> C d of ©, where t £ U and 

p(*o) = zo- 

3.1 Complex flow 

Given a path 7, with 7(0) = to G C and 7(1) = ii £ C. Let us assume that for every y £ N, 
where N is a neighborhood of xo £ C d , there exists a local solution tp(t ,y ) to the (CIVP) 

i(t) = f(t,x{t)), x{t ) = y eC d . 

For every such solution (p^ ,yo) ^ us assume, that the analytic continuation of <P(t ,y ) along 
7, denoted by <f(t ,y )i exists. In this situation, we have a map 

$}:N^C d , yo»$}yo:=<P { t ,y )(ti). 

Note As <3?J maps the initial value yo at to to the corresponding "solution value" ^yo at 
t%, we call this map the complex flow (or the evolution induced by /). The complex flow can 
be interpreted as the path-dependent function, mapping an initial value yo to the value of the 
analytic continuation <P(t Qj y ) evaluated at t±. 

3.2 Time grids and discrete flow 

In review of the introductory example, our presented real and complex approach took both 
use of discretizations of a given path 7 : [0, 1] — > U C C. Since C is not equipped with an order 
relation, we adapt the real-valued concept of consecutive points in time to < • • < t n G R, by 
the use of indexed sets 

A:= [#,...,**] CC, 

the time grids with «a £ N time steps := t^ +1 — tj, j G {0, . . . , n& — 1}. For every time 
grid A, we denote 

ta := max It,- I 

ie{o,...,n A -i} 1 J 1 
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as the maximum step size of A. To increase readability, we will often omit the symbol A 
(that is, we will identify Tj = Tj, n = n A , etc.), if mathematical definiteness does not suffer. 
In the example of Section 2 our attempt was to construct a corresponding grid function 

x A ■ A -> C d , 1 1 ► x A (t), 

with 

x A {t) ps <^(i), Vt G A, 

where ip : U — > C d is a solution of ([T|) and A C [7. To do so, we have started by setting 
x A (to) := x o- Next, we have used Euler's idea of small "tangential" update steps to re- 
cursively compute the missing x A (tj)'s. To generalize this idea, we just have to replace the 
method for computing a new value x A (tj + i) from an already computed point x A (tj). For this 
purpose, we introduce a function 

* : V -»■ C d , (t,s,x) ' — * 

where V C C x C x C d is an adequate set. The variables s, t and x play the role of current 
time, next time and current value, respectively. At this point, one naturally assumes that 
ty t,s x is defined for every choice of (t, s, x) G C x C x C d , with (s,x) G ilt and \t — s\ being 
sufficiently small. Under these assumptions VP is often denoted as the discrete evolution in 
analogy to the complex flow. To put things together, we define 

XA(t j+ i) := ^ +1 ^x A (tj), Vj e {0, . . . , n A - 1}. 

For latter purposes, let 

^ A x := x A (t n ), 

in analogy to $ 7 xo- This is the terminal point reached when the discrete evolution is 
developed along the time grid A. 

3.3 Errors 

Given a time grid A := [to, • • • ,t nA ] C C, let x A be a grid function generated by a one-step 
method applied to ([2]). As we have already seen by the example in the previous section, the 
approximation process is expected to cause grid errors Ej := s A (tj), where j G {0, . . . ,%} 
and 

£ A : A -> C d , tj i ^ ^'* x - a?A(*i)- 

At this point, we define for all < i < j < n A , := ; where 7^ represents the 

traverse sequentially visiting ti, . . . ,tj. 

The grid errors are generated by the local inaccuracy - the consistence error - of the discrete 
evolution denoted bj@ 

e{t,x,r) := $* +T ''x - ¥ +T > l x. 

Moreover, these local errors interact with the sensitivity of $^' io to perturbations of x. In 
general they can be damped, amplified or fortunately be extinguished throughout the one-step 
recursion. A special variant of the latter case will be analyzed in what follows. 

4 As ^ is a discrete evolution, e(t,x,r) is defined for all (t,x) G fi/, if r becomes sufficiently small. 



7 



4 Theorem linear case 



In this section we will reveal the secrets of the introductory example above in the more general 
context given by the (CIVP) 

x(t) = Ax, x(t ) = x , (3) 

where A G C dxd , d G N and (to, xq) G C x C d . Motivated by several numerical case studies, 
the effort to construct the corresponding complex flow <l?*'* xo for an arbitrary choice of t G C 
up to a desired accuracy by using a Runge-Kutta method (RKM) of order p G N, seems 
to be heavily dependent on the detailed choice of 7 - the path of integration. Pay attention 
to the fact, that there is no path-dependence from the analytical point of view, as already 
mentioned. 



For every n G N, let 

A? 



7 (0), 7 ( -V-.^f^) ,7(i; 



n I \ n 



In the context of ([3]), the convergence theory of one-step methods assures, that the family of 
grid functions (x&~>) cN corresponding to the used (RKM) converges towards ip(t) := <fr t,to xo 
with order p for an arbitrary choice of 7, where 7(0) = to an d 7(1) = t. In the following we 
will show, that for every (RKM) of order p, there exists a path of integration 7* , such that 

< :=e A7 .(t) = 



Roughly speaking, the latter equation states, that there exists a path of integration 7*, with 
7*(0) = to an d 7*(1) = such that the error at the ending point is of order p + 1 - a 
superconvergence effectS Furthermore, 7* can be chosen in such a way, that |7*| < - — 
which means, that 7* is also of practical interest. 



4.1 Effects of complex conjugation 

Before we start to proof the main theorem of this section, sketched above, we take care of the 
second observation^ mentioned by our introductory example: 

For a certain class of complex time grids, one can ensure that the computed value at the 
ending point of the time grid is real- valued. Let us therefore have a look at the detailed 
structure of the discrete evolution ^t t + T ' t x of an arbitrary s-stage (eRKM) applied to problem 
([3]). If 2ljj, bi G C, i, j G {1, . . . , s}, are the coefficients of the (eRKM), it follows that 

q/t+^x = x + hh, 
i=i 

where 

i-l 

i=i 

5 Note that the maximum grid error has not to be of order p + 1. 
6 Compare page [5] 
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is the i-th stage of the used (eRKM), i G {1, . . . , s}. In fact, no matter what the coefficients 
of the (eRKM) are, if applied to the linear ODE ([3]), then the evolution of one time step 
r G C from an initial point x G can be expressed as a simple matrix multiplication Ms. 
Furthermore the matrix M can be expressed as a (matrix) polynomial in tA. The precise 
statement is captured by the following lemma. 

Lemma 4.1. Let s G N and ^ t+T ^x be the discrete evolution of an arbitrary s-stage (eRKM) 
given by the coefficients 21 G C sxs and b, c G C s applied to problem Then it holds that 

yt+T,t x = j p( r ^) X) vx G C d , 

where P is a polynomial (applied to the matrix tA) of degree s, with coefficients 

Po,--- ,p s G Q[2l 2 ,i, 213,1,213,2, • • • ,2l S) i, . . . ,2l S)S _i, bi, . . . , b s }. 

Proof. Let us first have a look at the stages fcj. By induction over i G {1, . . . , s}, we will show 
that 

ki = Pi{rA)Ax, 

where Pi is a polynomial of degree i — 1 with coefficients in Q[2l2,i, • • • , Sl^i, . . . , Sl^j-i]: Con- 
sidering the first stage of the (eRKM) we get 

ki = Ax = Px(tA)Ax, 

where Pi := 1. Next, let i G {2, . . . , s}. By using the induction hypothesis, it follows that 

i=i 

i-l 

= Ax + t^ %,jAPj(TA)Ax 

3=1 

As one can easily see, the expression in the bracket above is induced by a polynomial Pi 
satisfying the desired conditions. In conclusion, 



yt+^x = X + T ^2 bih = ( I + biPi(TA)rA 



x 



i=l \ i=l 



yields the desired claim. At this point we would like to make the reader aware of the fact, that 
there is no problem of commutativity in the matrix polynomials above, because of A n A m = 
A m A n , for all n,m G N . □ 
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Remark The polynomial P of the previous lemma is also known as the stability function of 
the underlying (eRKM). 

Definition 4.1 (symmetric time grid). Let A := [t^, . . . , i^J be a time grid. A is called 
symmetric, if there exists a permutation ir G S nA with ir 2 = id, such that 

'j '7r(j+l)-l' 

for all j E {0, . . . , «a — !}• 

Theorem 4.1. iei A« 6e a symmetric time grid with t^ n = to G K and t A ™ G K, where 
n := 77-a- /f x n := x^j {^n) ^ s constructed by using an s-stage (eRKM) with corresponding 
coefficients 21 G R sxs and b,c£K s applied to problem |2j) a/ong A^ ; w;/iere ^4 6 M dxd ; i/ien 

Xo e R d ^ x n G M d . 

Proof. Let Tj := rj^' 1 , for j G {0, . . . , n — 1}. By using Lemma l4~T| it holds that 

n— 1 

j=0 
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where P is a polynomial of degree s with real-valued coefficients. Therefore it is sufficient to 
show, that the product is in R sxs . First we observe, that P(tA) = P(tA) and P(tA)P{o~A) = 
P(<jA)P(tA), for all <j,t G C. Thus, 

71-1 

J=0 



where 7r G 5 n is the corresponding permutation to the symmetric time grid An, with 7r 2 = id. 
Pay attention to the fact, that for every M G C sxs , it holds that MM = ~MM => M~M G 
W xs . □ 

As a result, Theorem 14. 1 1 i ustifies the second observation of the introductory example in the 
more general context of a (CIVP) given by an explicit linear autonomous system of first order 
differential equations. There the time grid was chosen to be symmetric along a half circle of 
the complex plane. Theorem 14.11 states that if we start from a real point the terminal point 
will be real as well. 



n p & a) ■ n pfaA)pfa u+iyi A) 



je{0,...,n-l} 
j+l=7r(j'+l) 



jS{0,...,n-l} 
j+l<7r(j+l) 



n pfaV' n pfaA)pfaA\ , 

ie{o,.. n-i} i i { 1 ;""™7 1 \ } * ^ ' 



4.2 Main theorem linear case 

This section is dedicated to the third observation (accuracy gain at the ending point) of the 
introductory example. Let A be a time grid with starting point to £ C and ending point t G C. 
The main idea is, to ascribe the leading term of the error e n = £a(£) a t the ending point to 
the step sizes Tj. This is achieved by a direct (but hard-core) expansion of $''*°:eo — ^ a xq. 
By doing so, we derive a condition for the step sizes which enables us to increase the 
order of convergence at the ending point t. 

Theorem 4.2 (main theorem linear case). Let 7 be a curve, with 7(0) = to £ C and 7(1) = 
t G C. For every n G N, we define 8 n := t^j and 

Tj (n) := t£» - tf, j G {0, . . . , n - 1}. 

Furthermore let C G M + be a constant, such that 5 n < C^fL, where £(7) := Jq \ j{t)\ dt is 
the length 0/7. 

If e n denotes the grid error at the ending point t G C corresponding to a p-stage ( eRKM) 
of order p G N applied to problem (0) along A n , then 

e n 1 ^ +i A p+1 e^~ to)A 
lim = lim Tp V(^(n)) p — — — — x . (4) 

Proof. To increase readability we define for all n G N, 

V»i(n):=2^ — ^ — and w » : = 2^ — fci — ' 

fc=0 ' fc=p+l 
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where j € {0, . . . ,n — 1}. Furthermore let E := e'^ 7 )"" 4 " and x n := ^^xq. To begin with 
our reasoning, we have a look at 

En _ ^ to X -X n _ J_ 

SP - sP - XP 

"n <Jn "n 



1 

w 

1 



81 

n-l 



n-l 



n 



.3=0 

71-1 
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3=0 3=0 



X . 



By expanding the first product, we get 



£n _ J_ 

On On 



n— 1 n— 1 



^o^(n) JJ ^(n) + p(n) 



i=o ;=o 

^3 



^0, 



where 



(5) 
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nP 



nP \t-t \ p nP-E' 



which implies lim n _ >00 ^fe^xo = 0. In the next steps we show, that we can omit the missing 
ipj(n) (products) and the higher order terms of the Co>j(n)'s in ([5]) as n tends to infinity. By 
the use of 

n— 1 n— 1 n— 1 

} ] Uj (n) • x n - } ] uj (n) ipi (n) • x 
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it holds that 



_. n— 1 

lim = lim ^ Vwj(n) • 

n— >oo o£, n— >oo n„ i — * 



lim 

n— >oo 



i=0 
n— 1 



1 - 1 (^(n)^ 1 1 ^ ^ (r^n) 

%U (p+i)i fci 



0^) ' 



j=0 fc=p+2 



Taking into account that 



. n— 1 oo 



(t»A)' 



.„ — — /■•! 



n— 1 oo 



" |t-t |f E £ 
1 U| j=0k=p+2 

j 71— 1 

" | t - t0 |*.„2|> 



nP(CL( 7 ) P|| 
fc! • n k 

E 1 



|t-t | P n 

and lim n ^oo x n = &' to xo = e^~ to ^ A xo, shows finally the desired claim. 



□ 



Note (s-stage (eRKM)) What happens if the number of stages of the (eRKM) is larger 
than its order of convergence? In the case of an s-stage (eRKM) of order p < s £ N (notice 
that p > s is not possible in the explicit case) one can derive a theorem analogous to Theorem 
14.21 We have omitted this case for a better readability of the preceding proof. The only 
difference to the results above is given by a constant C G C sxs (the p + 1 coefficient of the 
RUNGE-KuTTA-polynomial P(tA)), such that 



n-l 



lim —p = lim 

n— >oo On n—>oo 



* ( r j( n )) p+i • e 



P+ 1 . e (t-t )A 



3=0 



Jp + Ty. 



c 



XQ. 



The corresponding proof is achieved by some additional estimations of higher order terms of 
Ijr, analogously to the ones seen above. 



Note (implicit (RKM)s) Analogously to the proof above, one can show, that for an s-stage 
implicit (RKM) method, the structure of the error expansion (|4|) is also of the form 



^ n— 1 

lim xp ^2 ( r i( n )) P+ " C(t,t ,A)x , 

n-*oo n . =Q 



where C(t,t ,A) G C sxs . 
Proposition 4.1. Theorem 14.21 yields 
n-l 

Yl ( T i( n )) P+1 = °> Vn > iV G N 

3=0 



=> lim -p 

n— >oo o n 



lim 

n— >oo 



f) 1 



0. 



In the context of the premises of Theorem 14.21 this statement provides us a condition, which 
increases the order of a (RKM) of order p£N applied to ([3]) from p to p + 1 at the terminal 
point. 



13 



4.3 Superconvergent paths 

With Remark 14.11 in mind, our next goal is obvious: we want to construct complex detours 
that provide us with super convergence for linear ODEs. We consider an s-stage (RKM) of 
convergence order p G N and use it as a method to compute an approximation of the complex 
flow $*''°xo corresponding to problem ([3]), where t G C. To achieve superconvergence we are 
interested in time grid families (A = [t^, . . . , t^ A ] ) ngN , with 

1. t = to and t = t A A (in other words J2]=o~ l T t = t ~ *°)' 

2. Y^t® 1 ( r j A ) = (this provides superconvergence) and 

3. ta < where D G M + (prerequisite of Theorem 14. 2|) 

Fortunately the second criterion can be easily interpreted geometrically. This helps us to get 
a canonical family of time grids for every p G N: 

We use the fact that all n-th roots of unity sum up to zero. Setting ( n (j) := e 27 ™», j G Z, 
it holds for every a G C, that 

n-l 

J>Cn(j) =0- 
j=0 

Now we choose the time grid A in a way, that (r ? A ) p+1 , j G {0, . . . , ua — 1}> becomes such a 
collection of roots of unity. As illustrated by Figure O n consecutive n(p+ l)-th roots of unity 
Cn(p+l)(^)> • • • ) Cn(p+i)(k + ( n — I))? are transformed by complex exponentiation z z p+1 to 
the set of n-th roots of unity. 




6.(4) 

Figure 3: Transformation for n = 5, p = 2 and A; = 3. 



Lemma 4.2. Let n,p G N, G Z and a G C. T/ien it ZioWs t/tot 

fc+(n-l) 

E Kn( P+ l)(j)) P+1 = 
j=k 
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Proof. It holds that 



fc+(n-l) fe+(n-l) ,-. fc+(n-l) 

£ («Cn( P+ i)(j)) P+1 = a p+1 £ (e 2 -^l) P =^ +1 £ e 2 - 

j=k j=k j=k 

n-1 fc+(n-l) 

= cr > e « + cr > e « 

j=k j=n 

n-1 fc-1 

j=k 1=0 
n-1 

j=0 



□ 



In conclusion, one way to build a superconvergent time grid A = [t^, . . . ] is to take 
time steps rj^ represented by n consecutive n(p+l)-th roots of unity, which have been suitable 
globally scaled and rotated by a factor a £ C, such that X^o" 1 T t = * ~~ (compare Figure 
HI). A proper choice for this factor is 



t-t 

a :- 



■c-^k+tn— 1) * 
L,=fc Cn(p+l)U) 




Figure 4: Counterparts for n = 5, p = 2 and fc = 3 (no global scaling and rotation). 



To achieve a superconvergent family of time grids, it is easier to focus on the time steps 
Tj^ as done above. For coding and proof, we now change from the time steps to the time grid 
elements t^, j 6 {0, . . . , n^}. 

As one can see by example in FigureHl a superconvergent choice of the time steps 7"j\ induces 
that the time grid elements tj are located on a circle segment, depending on the order p E N 
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of the chosen method, connecting to,t G C. Figure H] shows the situation for p = 2. There, 
the corresponding time grid elements lie on a third segment of a circle. Compare this to the 
half circle segment of the introductory example, where the explicit Euler method has been 
used (p = 1). 

In conclusion, if the order of the method is p G N, locating the time grid elements 
equidistantly on a (p + l)-th circle segment, connecting to,t G C, provides the desired super- 
convergence effect. The following theorem formalizes this idea. 

Theorem 4.3. Let 7* := 7^ t , where 



to-t 



, 1— 2x 

e v p+i 



cos 



7T 



P+l 



+ 



t + t 



(6) 



Furthermore, let n G N and x A7 * 6e £/ie grid function generated by a s-stage (RKM) of order 
p G N applied to T/ien 

t = 7 *(0) and i = 7*(l). 

Ftirt/iermore, 

$*>*°x -x A7 * (t) 
lim -=-£2 = 0. 

Proof. By the use of Euler's formula, a simple calculation shows that to = 7*(0) and t = 
7*(1). Moreover 



tp-t 

2i sin -£r 

p+i 



g ' p+i — g p+i 



so it is not hard to see, that 

.a; 



const 



t ,t,p 



g n(p+l) _ \ 



VjG {0,...,n-l}, 



VjG{0,...,n-l}, 



which implies^] 5 n G 0(n _1 ) and additionally the equidistance of the time steps T^ n . Finally, 
it is sufficient to show, that 



n-l / x p+l n-l 
A"' ■' 



3=0 



E(^-'f*V +1 =o. 

3=0 



It holds that 



n-l 



3=0 



p+i 



•+^> 



n— 1 

j g p+l n — g p+l n =0 

i=o 

n— 1 1 -1 

E/ 27Ti j 27TI 1 2tH j \ P+ + 

j g p+l n p+l n — g p+l 11 =0 

j=0 

n-l 



_2zLil \ P+ 1 

j r /■! J » - J j y 

3=0 



0. 



=:(*) 



7 Compare to criterion 3 on page 1141 
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Since (*) is the sum of the n-th roots of unity, the desired claim follows by using the expansion 
given by Theorem 14.21 □ 

The theorem above yields immediately the following 

Corollary 4.1 (superconvergence of (RKM)s). Given a Runge-Kutta method of order 
p £N applied to {3p along 7*. The approximation at the terminal point is of order p + 1. 

Note As complex conjugation do not affect the argumentation so far, the previous corollary 
holds also for 7* := 7*, where 7*(i) := 7*(i) for all t G [0, 1]. 

5 Nonlinear case 

As we have seen in the last section, Runge-Kutta methods are super convergent along certain 
paths of integration, if the right-hand side of the corresponding (CIVP) is linear. Having this 
fact in mind, it is natural to ask whether this feature transfers to the case of an arbitrary 
right-hand side /. Again our study was triggered by several numerical experiments. These 
experiments suggested an adapted version of the above superconvergence statement, that also 
holds locally in the nonlinear case. The main idea is to study the order of convergence not 
for nA — ► 00, as in the linear case, but for t — * to. Let now «a = k be a fixed number of time 
steps, k G N, and let h := t — to- 

We will see, that according to the integration along a time grid A induced by 7f 0ji , it holds 
that 

e n eO(\h\ p+1 ), 

in contrast to the linear case, where e n G O ^~^tt^ • R° u ghly speaking, the superconvergence 

effect is independent of \t — to | in the linear case. 

Figure [5] represents the situation in which a numerical integration along the discretization 
A is done by two Runge-Kutta iterations (k = 2). If we regard this path as a scaled and 
rotated copy of a normalized^! circle segment, we can express the two time steps of A as 
follows. 

ti — to = <J\h and t% — tf~ = 02/1, 
where o~\ and o<i G C are the two time steps of the grid A 2 0,1 . 

The general case, k > 2, behaves analogously. For I £ {1, ... , k}, let 07 be the Z-th time 

step of the normalized grid A fc 01 . Due to this fact, the corresponding scaled time steps can 
be expressed by aih, where h := t — to- 

5.1 Composition methods 

Now, we draw a connection between compositions methods and our idea of numerical integra- 
tion along complex paths. Let ^> T := xJ/°+ r >° be the underlying autonomous^ discrete evolution 

8 This means a (p+ l)-th circle segment 7^ connecting and 1, p 6 N. 
9 This is w.l.o.g. no restriction. 
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of the basic (RKM). If we introduce the fc-term composition method 

T h := y akh o ...o ^ aih , (7) 

p 

the process of constructing x k := xa (t£) from xq := x& (t^), where A := A^ '' 1 , is given by 

x k = T h x . 

By the super convergence conditions of Section I3~3| our a\,...,ak satisfy the two equations 

ax + . . . + a k = 1 and a{ +l + ... + a{ +l = 0. (8) 

These conditions are exactly the criteria for a classical theorem on composition methods. 
From [HI p. 39, Theorem 4.1] it follows, that T h is a one-step method of convergence order 
at least p + 1. 

Remark The order conditions ([8]) can be found in 0. There, some specific composition 
methods have been introduced by explicitly solving the above order conditions over the reals. 
[9] and [13] also stated methods, by solving these equations over the complex numbers. In 
contrast to them, our starting point has not been ([8]). The adaption of super convergent paths 
Jt t from the linear theorem, solves implicitly the sufficient order conditions. 

5.2 Iterations 

The method given by (|7|) is called a /c-term composition method. If the underlying discrete 
evolution ^ has got order p G N, it follows that T is a one-step method of order at least p+1. 
As concatenation of (RKM)s is a (RKM) again, there is no problem to regard T as a new 
basic (RKM). T can again be used in the same way to construction a composition method - 
this time of order at least p + 2. Applying this process iteratively, it is possible to generate 
methods of arbitrary order of convergence. 

We start with a (RKM) ^ of order p. Iteration, leads to methods T(? (parameter k omitted 
in this notation), which are recursively defined by 

:= ¥ h , 

x r+l ,— L r ' 
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where A r := A^ '* and r > 0. T r h is a method of order p(r). Here 

p(0) :=p, 
p(r + 1) := p(r) + 5, 

where g £ N - the gain of order of convergence - depends on ty. For example, if \& is the 
discrete evolution corresponding to a (eRKM), g = 1. If ^ is the discrete evolution of a 
symmetric method, g = 2. 

With our geometric picture in mind, we now can reduce the method T r to the basic method 
Vf, applied along a certain path. The previous recursion for T r implies, that T^xq = ^ r xo, 
where T is a recursively defined time grid depending on the parameters r (depth of recursion), 
k (number of time steps), p (order of and h (global step size). We exemplify T for several 
choices of r in Figure El As one can see, T becomes more and more a fractal-like structure. 
Furthermore the time steps of T comply with the coefficients a r j, j € {1, 2} in [SJ p. 5, "two 
term composition"] (k > 2 analogously). 



Remark As mentioned above our approach leads to the same class of methods, already 
discovered in [9] and |13j . In contrast to the construction of /c-term composition methods by 
solving the system of order conditions ([8]), we explicitly get these methods by using the basic 
(RKM) along special complex time grids. The iteration done in [9] and [13], to construct 
higher order methods based on a certain (RKM), transfers (in our geometric point of view) 
to an application of this (RKM) along the recursively defined time grid T. 



5.3 A complex orbit 

As a benchmark and illustration of the ideas we just introduced, we consider a classical 
example of celestial mechanics - the restricted three body problem. The corresponding system 
of differential equations 

X\ = X\ + 2X2 — A — = — /i 



( (xi+/x) 2 +x 2)3 ^( (xi _ A)2+x 2)3 

(9) 



„. „ X2 X 2 

X2 = X2 — 2X\ — fl — = — /i- 



^/((X 1+/ X)2+X|) 3 ^((X 1 -A) 2 +X|) : 



is motivated by the motion of a satellite with respect to the gravitation potential induced by 
the moon and the earth. Thereby p, := 0.012277471 is the ratio of the moon corresponding to 
the mass of the total system and p, := 1 — p. Furthermore, the mass of the satellite could be 
neglected and (xi(t), X2{t)) represents the vector of the satellite's coordinates (the motion 
stays in a plane) in terms of a "rotating" coordinate system, whose origin represents the 
gravitational center of the moon and the earth and in which both celestial bodies stay on 
predefined points on the xi-axis. 

Due to the american mathematician R. Arenstorf0, we know, that there exists a periodic 



10 Compare pQ. At this point we want to mention, that Arenstorf showed the existence of such orbits by 
analytic continuation of parameterized two-body problems, for which the exact solutions (conic sections) 
are explicitly known. 
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solution to the corresponding initial value problem given by the system and a suitable 
chosen initial value. As illustrated by Figure [7J the error at the terminal point (integration 
from to = to T := 17.065216560157960, where T is the period of the corresponding exact 
orbit) by the use of the explicit Euler method (discrete evolution denoted by \£) is more 
than 55 times larger than with the composition method Ti (k = 2). Here, both approaches 
have been calculated by using equidistant time grids with the same number of time steps. 

6 Outlook 

We have seen, that composition methods with complex coefficients is nothing else than nu- 
merical integration along special complex time grids using a basic (RKM). As mentioned in 
the introduction, numerical integration along complex time grids involves several problems in 
general. 

First, one has to handle monodromy effects. This means, if one uses the discussed complex 
methods, one has to be aware of the fact, that branch points close to the real axis can 
be circulated by the complex path of integration, which causes a switch of the calculated 
solution's branch. 

Moreover, one has to develop a computing framework, that enables the integrator to mimic 
the concept of global analytic functions. 

That all this effort is a worthwhile endeavor can be seen as follows. On the hand, the 
algorithmic effort of numerically solving an initial value problem by taking a complex detour, 
can be reduced in a significant way (compare [1]). Thereby the idea is to stay away from 
singularities. On the other hand, the extended viewpoint given by studying differential equa- 
tions over the complex numbers enriches the structural understanding of the solution. For 
example, it might be possible and of practical interest to solve boundary value problems over 
the complex numbers, in situations where a real-valued approach would be ill-posed. 
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Figure 6: Time grid r for r = 1, 2, 3, 4, 5 and 11 (g = h = 1 and p = k = 2). 



Figure 7: The "exact" Arenstorf orbit (period T = 17.065216560157960, corresponding initial value 
(xi(0),x 2 (0),ii(0),i 2 (0)) T = (0.9940, 0.0, 0.0, -2.001585106379080) T ) given by 100000 real-valued equidistant 
steps of Dopri5 (5-th order (eRKM), solid), approximation by 100000 equidistant real- valued steps of the explicit Euler 
method ^ (dotted) and real part of the approximation by 50000 (2 "micro" steps for each "macro" time step, this means 
k = 2) time steps of the complex explicit Euler method (dashed), induced by Ti. 



